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ABSTRACT 

An  investigation  into  the  limitations  on  the  enhanced  field-emitted  current  density 
in  a  fast-pulsed  (rise-time  =  ns),  high  voltage  (>  106  V),  1-inch  vacuum  diode  was 
conducted  using  a  computer  simulation  based  on  the  Fowler-Nordheim  equation. 

Oscillations  in  the  emitted  current  density  (due  to  the  change  in  the  amount  of 
space  charge  within  the  gap)  were  found  to  quickly  decay  into  a  final  steady-state  for 
the  voltages  applied.  Steady-state  values  for  a  wide  variety  of  work  functions,  electric 
field  enhancement  factors  (based  on  the  theory  that  "whiskers"  on  the  cathode  surface 
experience  varying  degrees  of  enhancement),  and  applied  potentials  were  compared 
to  two  benchmarks:  the  amount  of  current  density  required  to  explode  a  whisker  in 
<  10  ns  by  joule  heating  (JE  =  109  A/cm2);  and  the  Child-Langmuir  (C-L)  space- 
charge-limited  current  density. 

Steady-state  values  were  found  to  be  less  than  JE.  One  model  of  the  formation 
process  of  a  plasma  at  the  cathode  surface  requires  that  JE  be  met  or  exceeded  by 
the  steady-state  value.  Thus,  such  a  model  is  not  supported  by  this  project's  findings. 

The  C-L  limit  is  based  on  a  thermionic-type  emission  process.  As  only  pure  field 
emission  (i.e.,  no  thermionic  emission  included)  was  considered,  the  steady-state 
values  were,  in  all  conclusive  cases,  less  than  the  corresponding  C-L  limited  values. 
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I.   INTRODUCTION 

The  breakdown  of  the  surface  of  the  cathode  in  a  pulsed  vacuum  diode  has  been 
the  subject  of  much  investigation  in  recent  years.  The  existence  of  crater-like  cathode 
spots  after  breakdown  can  be  explained  by  unipolar  arcing,  but  the  phenomena 
occurring  prior  to  the  arcing  are  still  debated. 

One  model  contends  that  a  metallic  "whisker"  on  the  surface  of  the  cathode 
provides  a  source  for  the  enhanced  emission  of  electrons.  After  a  period  of  time  (less 
than  10  ns  for  a  fast-pulsed  diode),  the  whisker  explodes  (due  to  the  rapid  build  up  of 
energy  produced  by  the  emission  process)  into  a  dense  plasma  which  in  turn  drives 
the  unipolar  arc. 

This  project  will  attempt  to  show  that  field  emission  alone  (i.e.,  not  including 
thermionic  emission)  from  a  whisker  cannot  reach  high  enough  values  to  trigger  the 
explosion  mechanism.  A  computer  simulation  utilizing  the  Fowler-Nordheim  equation 
for  field  emission  will  be  used.  The  results  will  also  be  compared  to  the  well-known 
Child-Langmuir  space-charge-limited  current  density  for  a  vacuum  diode. 

A  look  at  the  dimensionless  variables  used  in  the  program  and  the  program  itself 
are  included  in  the  appendices. 


II.   BREAKDOWN  AND  EMISSION  THEORY 

A.  BACKGROUND 

For  over  six  decades  scientists  have  been  investigating  pulsed- voltage  vacuum 
diode  breakdown.  With  the  advent  of  newer  technologies,  shorter  and  shorter 
voltage-pulse  rise  times  have  become  possible,  with  values  ranging  from  10-6  seconds 
in  1928  to  10~9  seconds  today  [Ref.  1].  Corresponding  to  today's  very  short  rise  times 
has  been  the  discovery  that  breakdown  even  occurs  on  the  nanosecond  time  scale. 
The  effects  of  the  breakdown  process  on  the  surface  of  the  cathode  have  been  well 
documented  [Refs.  1-3],  and  a  theory  of  the  formation  of  cathode  spots  (Fig.  1)  via 
unipolar  arcing1  has  recently  been  proposed  [Refs.  2,3].  In  dispute,  however,  is  the 
process  by  which  the  plasma  driving  the  unipolar  arc  is  formed. 

Cathode  spots  are  characterized  by  the  formation  of  a  dense  quasi-neutral 
plasma  near  the  cathode  surface.  There  are  two  basic  schools  of  thought  concerning 
the  primary  mechanism  for  breakdown  and  the  source  of  the  plasma,  only  one  of 
which  will  be  attended  to  here.  Anode-initiated  breakdown  would  occur  when  ions 
released  from  the  anode  (due  to  bombardment  of  the  anode  by  electrons  emitted  from 
the  cathode)  interact  with  the  cathode  to  form  the  plasma.  Here  fast-pulse 
(  =  nanoseconds)  breakdown  shall  be  investigated,  in  which  the  anode-emitted  ions 
would  not  have  time  to  reach  the  cathode  (for  the  1  in  =  2.54  cm  diode  gap  of  interest) 
before  breakdown  occurs;  anode-initiated  breakdown  shall  be  discussed  no  further.2 


1  The  concept  of  an  arc  which  originates  and  terminates  on  the  same  surface  (i.e.,  "unipolar") 
was  first  expounded  by  Robson  and  Thonemann  in  1958,  and  was  later  expanded  by  Schwirzke 
in  the  1980's[Refs.  2,4]. 

2  For  a  more  thorough  treatment  of  anode-initiated  breakdown,  see  [Ref.  1]. 
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Figure.  1.  Crater-like  cathode  spots  formed  by  unipolar  arcing. 


In  the  cathode-initiated  breakdown  school,  it  is  proposed  that  all  of  the  ions  and 
electrons  which  make  up  the  plasma  are  produced  by  processes  which  are  initiated  by 
events  occurring  only  at  the  cathode  [Ret.  1].  Within  this  school,  however,  there  are 
many  camps,  each  with  its  own  views  on  what  processes  do  indeed  occur,  and  which 
events  really  take  place  to  initiate  them.  One  of  the  more  reknown  camps  is  led  by 
Schwirzke,  who  recently  proposed  that  neutrals  released  from  the  surface  of  the 
cathode  are  ionized  by  the  emitted  electrons,  thus  forming  a  plasma  near  the  cathode 
surface  [Refs.  2,3].  Another  is  led  by  Dyke  (et  al)3  who  proposed  that  joule  heating  of 
an  electron  emission  site  causes  the  site  to  vaporize,  and  in  the  process  ions  are 
created  which  allow  the  plasma  to  form  [Ref.  1].  Here  shall  be  examined  a  third  camp, 
championed  by  Mesyats  [Ref.  1],  in  which  Dyke's  theory  is  taken  to  its  most  extreme 
limit,  that  of  an  explosion  of  the  emission  site  into  a  plasma  directly.  The  goal  is  to 
show  that  an  elementary  version  of  a  Mesyats-type  "explosive"  model  cannot  support 
breakdown  on  a  time-frame  indicative  of  fast-pulsed  voltage  vacuum  diode  breakdown. 

B.  THE  MODEL  DIODE 

The  model's  anode  and  cathode  will  consist  of  two  polished  metal  plates 
separated  by  a  vacuum.  Despite  the  polishing  and  the  degassing  effect  of  the 
vacuum,  however,  some  small-scale  "irregularities"  will  nevertheless  exist  on  the 
surfaces  of  the  plates.  Weakly  bound  bits  of  dust,  oxides,  and  non-metallic  inclusions 
litter  the  cathode  surface,  providing  excellent  sites  for  the  enhanced  field  emission  of 
electrons  [Ref.  2].  These  so-called  "whiskers"  will  vary  in  size  and  shape,  with  the 
more  needle-like  (height-to-base  ratio  »  1)  metallic  protrusions  providing  the  most 


3  See  [Ref.  1]  for  a  list  of  references  on  Dyke's  experiments. 


enhancement  of  the   field   emission.     Flatter,   non-metallic   (e.g.,   dielectrics,   semi- 
conductors, etc.)  sites  can  also  provide  significant  enhancement  [Ref.  5]. 

C.  ELECTRON  EMISSION  FROM  THE  CATHODE 

There  are  two  mechanisms  for  the  electric-field-enhanced  emission  of  electrons 
from  a  metallic  surface:  thermionic  (high  temperature,  low  applied  field);  and  field  (low 
temperature,  high  applied  field)  [Ref.  4].  The  presence  of  the  electric  field  at  the 
surface  lowers  the  potential  barrier  at  the  plate-vacuum  interface,  thus  allowing 
electrons  in  the  conduction  band  an  opportunity  to  become  free  from  the  metal  (Fig.  2). 
In  thermionic  emission,  electrons  at  high  temperatures  gain  enough  kinetic  energy  to 
pass  over  the  barrier  and  into  the  gap.  In  field  emission,  the  electrons  tunnel  through 
the  barrier  instead.  In  general,  the  electric  field  at  the  emitter  surface  must  be  on  the 
order  of  107  V/cm  for  field  emission  to  take  place  [Ref.  6].  For  a  diode  gap  of  2.54  cm 
(=  1  in),  an  applied  potential  of  2.54x1 07  V  would  be  necessary  to  achieve  field 
emission,  but  as  noted  previously,  the  whisker  emission  sites  experience  an  enhanced 
applied  electric  field  (Fig.  3).  The  enhancement  factor,  m,  varies  from  1  to  several 
hundred,  with  needle-like  whiskers  having  the  largest  values  [Ref.  7].  For  electric 
fields  less  than  107  V/cm,  thermionic  emission  will  dominate  for  temperatures  above 
approximately  1500  K,  with  little  substantial  emission  at  lower  temperatures  [Ref.  6]. 

For  the  case  of  the  non-metallic  inclusions,  a  combination  of  field  and  thermionic 
emission  occurs.  Electrons  tunnel  from  the  metallic  cathode  material  into  the 
dielectric,  then  are  thermionically  emitted  into  the  gap  (see  Fig.  4)  [Ref.  5]. 
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Figure  2.   Lowering  of  potential  barrier  at  metal-vacuum  interface  [Ref.  7]. 


E-FIELD   LINES 


nTTy? 


~10_4cm 


////// 


Figure  3.   Enhancement  of  electric  field  at  whisker  tip  [Ref.  7]. 


VACUUM 


Bortiolly  emoroeee 
"on-mtlollic  impurity 


SUBSTRATE   /    /  / 
METAL  CATHODE 
/////// 


0»cr«stmq 


METAL 


Figure  4.  Combination  of  field  and  thermionic  emission  for  a  dielectric  [Ref.  5]. 

For  this  elementary  model,  only  field  emission  from  metallic  whiskers  shall  be 
considered,  with  the  following  assumptions: 


1.  The  time-scale  for  the  explosion  of  a  whisker  is  only  a  few 
(<  10)  nanoseconds; 

2.  The  applied  voltage  rises  instantaneously  to  a  constant  value  at  t  =  0; 

3.  The  gap  is  perfectly  evacuated;  no  positive  ions  are  present,  so  no  heating  of 
the  site  by  ion  bombardment  is  possible; 

4.  The  initial  temperature  of  the  cathode  (=300K)  is  too  low  for  any 
appreciable  amount  of  thermionic  emission  to  occur  at  first;  and 

5.  Any  increase  in  temperature  during  the  short  time-frame  prior  to  explosion 
shall  be  neglected  (and  thus  so  shall  thermionic  emission)  in  this  elementary 
model. 


Pure  field  emission  is  governed  by  the  Fowler-Nordheim  equation  [Ret.  6] 

F2       (  -6.83x1 07)f(V)V 

J(F)  =   (  1.54X10-6)  —  e  F     .  (1) 

Since  f(y)  (the  Nordheim  elliptic  function)  is  slowly  varying,  it  can  be  approximated  as 
a  constant4 

f(y)  =  0.94  , 

which  leads  to  (see  Figs.  5a,b) 

A3'2 

_       c2       ( -6.42x1 07)-S=- 

J(F)  =   (  1.54x10"*  J-^-e  F     ,  (2) 

where 


J  ■  emitted  current  density  (A /cm2  ); 

F  =  -mE  s  enhanced  electric  field5  at  the  whisker  tip  (V/cm  ); 

m  ■  electric  field  enhancement  factor;  and 

4>  =  work  function  of  the  cathode  metal  (eV). 

J  is  a  function  only  of  F  for  any  given  whisker;  any  negative  (i.e.,  electron)  space 
charge  which  exists  in  the  gap  will  reduce  the  emitted  current  density.  If  the  applied 
potential  remains  constant,  it  can  be  predicted  that  J  will  start  at  some  initial  value,  J0  , 
and  decrease  as  more  and  more  space  charge  enters  the  gap.  The  electrons,  once  in 
the  gap,  will  be  accelerated  across  to  the  anode,  eventually  disappearing  there.  When 
the  first  emitted  electrons  reach  the  anode,  the  electric  field  at  the  cathode  will 
increase  due  to  the  loss  of  space  charge.  The  increased  field  will  cause  an  increase 


4  For  a  more  thorough  treatment  of  f(y),  see  [Refs.  4,6]. 

5  Eq.  (2)   requires  that  F  be  defined  as  a  positive  quantity,  hence  the  negative  sign. 
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Figure  5a.   Log  plot  of  J  vs.  F  (Eq.  2). 
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Figure  5b.   Expanded  log  plot  of  J  vs.  F  for  various  work  functions  (Eq.  2). 
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in  J,  but  soon  the  increased  addition  of  space  charge  will  exceed  the  loss  at  the 
anode;  the  amount  of  space  charge  in  the  gap  will  then  rise,  the  electric  field  at  the 
cathode  will  decrease,  and  J  will  fall.  This  oscillatory  action  should  diminish  as  the 
emitted  current  density  at  the  cathode  comes  to  equal  that  going  into  the  anode.  In 
other  words,  the  current  density  becomes  constant  at  some  final  space-charge- 
determined  value,  JF. 

It  should  be  noted  that  in  reality,  the  applied  potential  is  not  constant  from  t  =  0. 
Rather,  it  ramps  up  from  zero  to  some  final  value.  The  ramping  has  the  effect  of 
causing  the  current  density  to  rise  in  a  similar  fashion;  oscillations  such  as  those 
described  above  are  not  observed. 

D.  NON-APPLICABILITY  OF  CHILD-LANGMUIR  LAW 

It  is  somewhat  important  to  note,  at  this  point,  that  JF  is  not  the  well-known 
Child-Langmuir  space-charge-limited  current  density  (A/cm2) 


Jcl  =    g-  Co 


2e 
M 


1/2       Q/0 
\/3/2 


(3) 


with 


£o  =  permitivity  of  free  space  =  8.85x1 0"12  C2/N-m2  ; 
e  =  I  electron  charge  I  =  1 .60x1 0~19  C  ; 
M  ■  electron  mass  =  9.1 1x1 0"31  kg  ; 
V  =  applied  voltage  (V)  ;  and 
d  a  diode  gap  length  (cm)  . 

JF  should  be,  in  fact,  somewhat  less  than  JCL  .  This  difference  arises  due  to  the  fact 
that  JCL  is  not  derived  for  the  case  of  a  pure  field  emission  mechanism.  Rather,  there 
must  be  an  ample  supply  of  electrons  (eg.,  in  the  case  of  electrons  emitted  from  a 
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heated  filament)  always  available  for  emission.  In  the  derivation  of  JCL,  one  of  the 
boundary  conditions  is  that  the  electric  field  at  the  cathode  go  to  zero  due  to  the 
presence  of  the  space  charge  (see  [Refs.  4,6,7]).  If  this  were  to  happen  to  a  pure  field 
emitter,  then  J  would  go  to  zero,  and  a  steady  state  would  not  develop.  Since  the 
field  emitter  cannot  emit  at  JCL  for  any  appreciable  length  of  time  without  shutting  itself 
off,  its  steady  state  emission  value  must  be  something  less: 

Jf  <  jcl  ( Pure  *'e|d  emission  )  .  (4) 

E.  EXPLOSIVE  MODEL  CURRENT  DENSITY  REQUIREMENT 

It  has  been  shown  by  Schwirzke  [Ref.  2]  that,  for  stainless  steel,  a  current  density 
of  JE  =  109  A/cm2  would  be  required  to  explode  a  whisker  on  the  cathode  surface  in  a 
few  (<10)  ns  by  joule  heating.  Thus,  a  whisker  would  need  to  emit  a  Javg  =  JE  for 
approximately  10  ns  before  exploding.  If  the  transient  oscillation  period  is  small 
compared  to  10  ns,  then  Javg  =  JF,  and 

JF  >  JE  (explosive  emission)  (5) 

is  required  for  the  explosive  emission  model.  This  project  will  attempt  to  show  that,  for 
a  wide  range  of  applied  potentials,  work  functions,  and  enhancement  factors,  Eq.  (5) 
cannot  be  satisfied,  thus  proving  that  field  emission  alone  is  not  capable  of  providing 
the  current  density  required  by  the  explosive  model. 
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III.  SIMULATION 

A.  BACKGROUND 

The  type  of  simulation  model  used  here  is,  as  [Ret.  8]  so  aptly  puts  it,  one  of  the 
"simplest"  possible,  viable  since  the  first  digital  computers  came  into  use  in  the  1950's. 
While  almost  primitive  in  comparison  to  the  latest  innovations  in  so-called  "particle-in- 
cell"  simulations  [Ref.  8],  the  model  used  here  nevertheless  is  time-tested  and  has 
worked  quite  well  for  the  rather  simple  events  under  investigation.  Although  the  basic 
techniques  are  over  three  and  a  half  decades  old,  the  actual  program  is  original  and  is 
included  (in  one  form)  in  Appendix  B  for  scrutiny.  The  program  is  written  in  the  "C" 
language  and  was  run  on  a  Sun  SPARCstation  2  incorporated  into  a  system  running 
SunOS™  4.1. 1.. 

B.  SIMULATION  THEORY 

In  general  terms,  the  simulation  works  as  follows.  Space  charge  is  simulated  not 
as  a  continuous  beam  of  electrons,  but  rather  as  a  series  of  infinitely  thin,  numerous 
charged  disks  of  varying  charge  density  o  [C/m2].  In  each  time  interval,  the  forces  on 
the  disks  within  the  gap  are  computed,  the  disks  are  moved,  and  a  new  disk  is  created 
at  the  cathode.  When  a  disk  reaches  the  anode,  it  is  no  longer  used  in  any  further 
calculations.  In  the  final  steady  state,  all  of  the  disks  will  have  the  same  charge 
density;  for  every  time  step  just  as  much  charge  is  lost  at  the  anode  as  is  created  at 
the  cathode. 
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1.  Electric  Field 

The  charge  density  a  is  calculated  by  multiplying  the  time  step  At  by  the 
current  density  J  (Eq.  2)6 

o  =  -J  At  .  (6) 

The  field-emitted  current  density  at  the  cathode  surface  depends  on  E^O), 
the  value  of  the  electric  field  at  the  cathode  due  to  all  of  the  space  charge  disks,  and 
Eapp,  the  component  of  the  electric  field  at  the  cathode  surface  due  to  the  applied 
potential  across  the  diode.  In  order  to  find  Esc(0),  one  starts  with  the  one-dimensional 
Poisson's  Equation  [Ref.  6] 

d  2  V(x)         q(nj-ne) 
dx^  £o 


which,  for  the  ion  density  n,  =  0  and  the  electron  charge  q  =  -e  becomes 

(8) 


d  2  V(x)         en(x) 


dx^  £o 

where 

V  ■  potential  (V)  ; 

x  =  position  (  =  0  at  cathode,  =  d  at  anode  )  (m)  ; 

e  =  I  electron  charge  I  =  1 .60x1 0"19  C  ; 

n  ■  electron  number  density  (m-3)  ;  and 

e0  ■  permitivity  of  free  space  =  8.85x1 0~12  C2/N-m2 


6  The  current  density  J  will  always  be  defined  as  positive,  hence  a  negative  sign  has  been  ad- 
ded on  the  right  side  of  Eq.  (6). 
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Using  dV(x)/dx  =  -E(x),  one  gets 

dE(x)    B    -en(x)  (g) 

dx  £o 

Integrating    both    sides    over   the   gap,    using    the    relation    Esc(d)  =  -Esc(0),    and 
rearranging  yields 

d 


E 


(0)  =  -l-Joen(x)dx  .  (10) 


2eo 


Now,  using  the  charged  disk  concept,  the  integral  becomes  a  summation 

Esc(0)  =  -=J-  ZenkAx   .  (11) 

"0   gap 

where        nk  ■  electron  number  density  [m-3]  in  disk  number  k.  Realizing        that 

e  n  Ax  =  -o  yields 

Esc(0)  =  ■=*-  Z  -<*k   •  (12) 

^MD    gap 

Finally,  using  Eq.  (6), 

Esc(O)  =  ^-  Z  JkAt  .  (13) 

"0   gap 

For  a  disk  located  at  a  point  x'  in  the  gap,  Esc(x')  can  be  determined  from 
Eqs.  (9)  and  (10)  by  breaking  the  integral  of  the  left-hand-side  into  two  parts  and 
subtracting: 

J  dEsc(x)  =  —I  -en(x)dx  (14a) 

becomes 

Esc(x')-Esc(0)  =  — 1-   Z  4At  ,  (14b) 

Eo  o-*x' 
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and 


Esc(d) 


d 


becomes 


J   dEsc(x)  =  -1  J   -en(x)dx  (15a) 


Esc(d)-Esc(x')  =  --1   2  JkAt  .  (15b) 

£0    x'-xJ 

Subtracting  Eq.  (15b)  from  Eq.  (14b)  and  rearranging  yields 

Esc(x')  =   ^-    Z   JkAt  -   -J-    £  JkAt   .  (16) 

^   x'-»d  «0   0-*x' 

Eq.  (16)  is  equivalent  to  saying  that  the  space  charge  "in  front"  (x'-kj)  will  contribute 
positively  to  the  electric  field  at  x',  whereas  the  space  charge  "behind"  (0-»x')  will 
contribute  negatively. 

If  the  applied  voltage  remains  constant,  so  will  the  applied  electric  field  Eapp. 
The  total  electric  field  at  x'  is  therefore7 

EtOO  =  -Eapp  +  77   IJkAl-TT   I  JkAt.  (17) 

d£C     X'-Mj  "0     0-»X' 

At  the  cathode,  the  enhancement  factor  must  be  included: 

F  =  -m  ET(0)  =  m  JEapp  -  ^-  £  JkAt|  .  (18) 

A  schematic  diagram  of  the  electric  field  components  in  the  diode  is  shown  in  Fig.  6. 


7  Like  F,  Ew  will  be  defined  as  a  positive  quantity,  so  the  appropriate  sign  change  must  be 
made. 
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+ 


^ — > 


Figure  6.  Schematic  diagram  of  electric  field  components  in  diode. 


17 


2.  Equation  Of  Motion 

The  nature  of  the  large  (>  106  V)  applied  potentials  used  in  the  simulation 
requires  that  relativistic  equations  of  motion  be  used.  For  any  electron  in  a  disk  at  a 
position  x,  the  relation  is  (based  on  mass  x  acceleration  =  charge  x  E-field) 


where 
ys  1  /  Vl  -  v2  /  c2  ; 
P  =  V/c; 
t  ■  time  (s)  ; 

e  =  I  electron  charge  I  =  1 .6x1 0"19  C  ; 
M  =  electron  mass  =  9.1 1x1 0"31  kg  ; 
c  =  speed  of  light  in  vacuum  =  3.0x1 08  m/s  ; 
E\-(x)  s  total  electric  field  at  position  x  (V/m)  ; 
x  a  position  (  =  0  at  cathode,  =  d  at  anode  )  (m)  ;  and 
V  ■  velocity  (m/s) . 

Substituting  into  Eq.  (19)  and  carrying  out  the  derivative  gives 


Taking  the  derivative  of  y  and  substituting  in  yields 


aV        -e 


With  a  little  more  substituting,  Eq.  (21)  can  be  solved  for  dv*/  dt: 


(ft         -e 


dt         My3 

Eq.  (22)  is  the  equation  of  motion  used  in  the  simulation 
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(19) 


^♦Tf-T?***-  <2°> 


(PV  +  Y)^-  =   -^t(x)   •  (21) 


£T(x)   .  (22) 


IV.   RESULTS 

A.  OBJECTIVE 

Again,  the  goal  of  the  simulations  is  to  show  that  the  value  of  Jp  for  a  given  set 
of  conditions  will  not  be  greater  than  the  109A/cm2  required  by  the  elementary 
explosive  emission  model  (Eq.  (5)). 

B.  METHODS  OF  OBTAINING  RESULTS 

For  lesser  values  (<  106  A /cm2)  of  J0  (the  value  of  the  current  density  at  at  the 
cathode  at  t  =  0;  i.e.,  with  the  potential  applied  but  no  space  charge  in  the  gap  yet) 
the  program  is  easily  capable  of  running  through  the  transient  period  and  into  the  final 
state  (as  in  Fig.  7a)  in  a  short  (=  computer  hours )  period  of  time.  But  for  larger 
values  of  J0,  the  run  time  goes  up  considerably  (=  computer  days ).  The  objective 
requires  only  that  JF  <  JE  is  proven;  it  is  possible  to  do  this  without  running  the 
program  through  to  the  final  state  by  using  three  simple  "tricks": 

1.  Attempt  to  simulate  the  steady-state  electric  field  at  the  cathode  by  artificially 
adding  in,  at  t  =  0,  an  estimate  of  the  steady-state  space  charge.  The 
artificial  charge  will  be  an  immovable  source  located  at  the  anode  so  that  it 
will  always  be  "in  front"  of  the  emitted  charges; 

2.  Remove  equal  parts  of  the  artificial  space  charge  during  each  time  step  of  the 
simulation,  until  it  is  all  removed;  and 

3.  Set  the  amount  of  artificial  space  charge  removed  each  time  step  equal  to  the 
amount  created  at  the  cathode  in  the  first  time  step. 

Although  the  resulting  plots  will  certainly  be  non-physical  (no  such  artificial  space 
charge  exists  in  reality),  the  graphical  output  can  nonetheless  be  definitively 
interpreted. 
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During  the  early  phase  of  the  removal  period,  the  graph  of  J  vs.  t  will  be  a 
horizontal  straight  line,  as  in  the  final  steady  state  (just  as  much  space  charge  is 
removed  each  period  as  is  put  in,  so  J  will  not  change)  (see  Figs.  7b,c,d).  If  the 
estimate  of  the  steady-state  space  charge  is  exactly  correct,  the  graph  of  J  vs.  t  will 
continue  to  be  straight,  even  after  all  of  the  artificial  space  charge  is  removed, 
because  in  this  case  the  simulation  will  have  started  equivalent  to  the  steady  state. 
Thus,  if  the  graph  remains  flat  after  all  the  artificial  space  charge  is  gone,  the  value  of 
JF  will  be  known  exactly  (see  Fig.  7b).  Guessing  the  exact  value  of  the  steady  state 
space  charge  is,  however,  highly  unlikely. 

If  the  estimate  of  the  steady-state  space  charge  is  too  high,  the  value  of  F  will  be 
smaller  than  in  the  actual  steady-state  case,  and  the  value  of  J  during  the  early  phase 
of  the  removal  period  will  be  lower  than  the  actual  JF.  In  this  case  some  of  the 
artificial  charge  will  still  exist  when  the  first  emitted  electrons  reach  the  anode.  At  this 
point,  the  amount  of  charge  lost  each  time  period  will  increase,  thus  increasing  F  and 
J  (see  Fig.  7c). 

More  useful,  however,  is  the  case  where  the  estimate  of  the  steady-state  space 
charge  is  too  low.  The  value  of  F  will  be  larger  than  in  the  actual  steady-state  case, 
and  J  will  be  larger  than  the  actual  JF.  In  this  case,  the  artificial  charge  will  be  used  up 
before  the  first  emitted  electrons  reach  the  anode.  Once  all  of  the  artificial  charge  is 
gone,  the  value  of  F  will  drop  as  more  charge  enters  the  gap  while  none  is  removed 
(see  Fig.  7d).  It  can  be  concluded  that  the  early  removal  phase  value  of  J  is  greater 
than  the  actual  JF.  An  upper  bound  on  JF  can  thus  be  found  by  this  "maximum  value 
method",  and  if  the  calculated  upper  bound  is  less  than  JE,  the  objective  will  have 
been  met. 
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Figure  7a.  Log  of  the  current  density  ratio  vs.  time.  No  artificial  space  charge  is 
added,  and  the  simulation  is  run  until  the  final  steady  state  is  reached 
(2  ns).  The  steady-state  ratio  is  J/J0  =  1.37x10"2. 
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Figure  7b.  Log  of  the  current  density  ratio  vs.  time.  Enough  artificial  space  charge 
was  added  to  simulate  the  steady  state  of  Fig.  7a  exactly  (12.3752 
dimensionless  units).  The  value  of  f  is  the  number  of  time  steps  required 
to  use  up  the  artificial  charge.  The  value  of  old  indicates  the  time  step  of 
creation  for  the  charge  disk  closest  to  the  anode.  If  f  >  0,  old  >  1 ,  and  the 
plot  is  straight,  then  the  line  is  at  the  equilibrium  value.  Note  the  run-time 
is  0.1  ns  (1/20  of  the  full-run  time-steps  in  Fig.  7a). 
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Figure  7c.  Log  of  current  density  ratio  vs.  time.  More  than  the  equilibrium  value  of 
the  space  charge  has  been  artificially  inserted  (14  dimensionless  units). 
In  this  case,  f  =  0,  old  >  1  indicates  that  the  artificial  space  charge  has  not 
been  used  up  before  the  first  disk  reached  the  anode;  the  plot  rises  at  the 
end  as  expected.  Again,  only  1/20  of  the  full-run  time-steps  is  required 
to  obtain  a  result  (in  this  case,  a  lower  bound  on  JF). 
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Figure  7d.  Log  of  current  density  ratio  vs.  time.  Less  than  the  equilibrium  value  of 
the  space  charge  has  been  artificially  inserted  (11  dimensionless  units). 
In  this  case,  f  >  0  and  old  >  1 ,  but  the  plot  drops  from  the  early-phase  flat 
portion.  Note  that  the  early  phase  value  (-1.588)  is  greater  than  the 
steady  state  value  (-1.863)  of  Fig.  7a.  An  upper  bound  on  the  current 
density  has  been  determined  in  only  1/20  of  the  full-run  time-steps. 
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C.   SIMULATION  RESULTS 

As  can  be  seen  in  the  last  column  of  Table  1 ,  under  nearly  all  tested  conditions, 
JF  <  Je  =  !09  A/ cm2.  In  fact,  in  most  cases,  J0  was  less  than  JE,  indicating  that  only 
for  high  applied  electric  fields  (>  3x1 06  V/cm)  and  large  enhancement  factors  (>100) 
would  there  even  be  a  possibility  for  the  final  value  of  the  current  density  to  exceed 
the  explosive  limit. 

For  those  cases  where  it  was  not  shown  that  JF  <  JE  (denoted  by  t  ).  the 
excessive  amount  of  run  time  ( =  computer  days )  to  lower  the  maximum  level  was  just 
too  prohibitive. 

It  was  predicted  earlier  that  JF  would  be  less  than  the  Child-Langmuir  space 
charge  limited  current  density  JCL  (Eq.  (4)).  The  program  was  run  long  enough  (when 
possible)  for  the  <J>  =  3.0  eV  case  to  determine  if  the  maximum  of  JF  is  indeed  less  than 
JCL.  The  results  are  summarized  in  Table  II.  Similar  results  are  expected  for 
<l>  =  4.0  and  5.0  eV  as  well,  since  the  values  of  J0  (and  logically  the  maximum  JF  also) 
are  lower  for  these  two  cases  than  for  the  <t>  =  3.0  eV  case. 
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TABLE  I.   SIMULATION  RESULTS t 


EaDD(V/cm) 

0(eV) 

m 

J0  (A/cm2) 

J^(A/cm2) 

3.94  x105 

5.0 

1 

=  0 

=  0 

3.94x105 

5.0 

10 

=  0 

=  0 

3.94  x105 

5.0 

100 

5.86 

4.93 

3.94  x105 

5.0 

300 

9.92x1 06 

£  7.84x1 06  * 

3.94  x105 

4.0 

1 

=  0 

=  0 

3.94  x105 

4.0 

10 

=  0 

=  0 

3.94x105 

4.0 

100 

1.30x103 

8.69x1  O1 

3.94  x105 

4.0 

300 

6.98x1 07 

<  6.97x1 07  * 

3.94  x105 

3.0 

1 

=  0 

=  0 

3.94  x105 

3.0 

10 

=  0 

=  0 

3.94  x105 

3.0 

100 

1.68x105 

2.17x102 

3.94  x105 

3.0 

300 

4.27x1 08 

<  5.00x10s  * 

3.94  x106 

5.0 

1 

=  0 

=  0 

3.94  x106 

5.0 

10 

5.86 

5.77 

3.94  x106 

5.0 

100 

7.73x1 09 

<6.81x107* 

3.94  x106 

5.0 

300 

2.34x1  O11 

<  4.24x1 07  * 

3.94  x106 

4.0 

1 

=  0 

=  0 

3.94  x106 

4.0 

10 

1.30x103 

4.84x1 02 

3.94  x106 

4.0 

100 

1.62x1010 

<  7.24x1 07  * 

3.94  x106 

4.0 

300 

3.48x1  O11 

<  7.39x1 07  * 

3.94  x106 

3.0 

1 

=  0 

=  0 

3.94  x106 

3.0 

10 

1.68x105 

2.30x1  O3 

3.94  x106 

3.0 

100 

3.42x1  O10 

<  8.00x106  * 

3.94  x106 

3.0 

300 

5.41x1011 

<9.00x106  * 

3.94  x107 

5.0 

1 

5.86 

5.85 

3.94  x107 

5.0 

10 

7.73x1 09 

<6.81x107  * 

3.94  x107 

5.0 

100 

3.98x1 012 

<  8.00x1 08  * 

3.94  x107 

5.0 

300 

4.05x1 013 

<  7.00x1 09*  * 

3.94  x107 

4.0 

1 

1.30x103 

1.07x103 

3.94  x107 

4.0 

10 

1.62x1010 

<  7.24x1 07  * 

3.94  x107 

4.0 

100 

5.25x1 012 

<  8.00x108  * 

3.94  x107 

4.0 

300 

5.15x1013 

<  7.00x1 09*  * 

3.94  x107 

3.0 

1 

1.68x105 

1.55x104 

3.94  x107 

3.0 

10 

3.42x1  O10 

<  8.00x1 06  * 

3.94  x107 

3.0 

100 

7.32x1 012 

<  8.01x10s  * 

3.94  x107 

3.0 

300 

6.97x1 013 

<  7.00x1 09*  * 

t  For  a  diode  with  plate  spacing  of  1  in  ■  2.54  cm. 
*  Indicates  maximum  value  method  result. 
t  Indicates  run-time-constrained  result. 
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TABLE  II. 
SIMULATION  RESULTS  VS.  CHILD-LANGMUIR  LAW* 


+  (eV) 

m 

Eaoo  ( V/cm ) 

JCL  (A/cm2) 

J^A/cm2) 

maximum  JF  (A/cm2) 

3.0 

1 

3.94x1 07 

3.62x1 05 

1.68x105 

1.55x104 

3.0 

10 

3.94x1 06 

1.14x104 

1.68x105 

2.30x1 03 

3.0 

10 

3.94x1 07 

3.62x1 05 

3.42x1  O10 

8.00x10s* 

3.0 

100 

3.94x1 05 

3.62x102 

1.68x1 0s 

2.17x102 

3.0 

100 

3.94x1 06 

1.14x104 

3.42x1  O10 

8.00x10s* 

3.0 

100 

3.94x1 07 

3.62x1 05 

7.32x1 012 

8.01x108* 

3.0 

300 

3.94x1 05 

3.62x1 02 

4.27x1 08 

5.00x105* 

3.0 

300 

3.94x1 06 

1.14x104 

5.41x1011 

9.00x10s* 

3.0 

300 

3.94x1 07 

3.62x105 

6.97x1013 

7.00x1 09* 

D.   OTHER  GRAPHICAL  OUTPUT 

The  plot  of  Fig.  7a  appears  to  agree  with  the  basic  theory  that  the  value  of  the 
emitted  current  density  will  oscillate  at  first,  then  settle  down  until  a  final  steady  state  is 
reached.  But  earlier  versions  of  the  program  encountered  problems.  In  an  effort  to 
determine  the  causes  of  the  difficulties,  the  program  was  altered  to  produce  plots  of 
other  variables.  Figs.  8-1 1  are  some  examples  run  for  the  same  conditions.  Although 
the  applied  electric  field  is  rather  high  in  these  plots,  the  results  are  nonetheless 
representative  of  the  output  for  the  applicable  range  of  electric  fields.  Short  computer 
running  time  was  the  primary  factor  in  the  choice  of  parameters  for  these  examples. 


t  For  a  diode  with  plate  spacing  of  1  in  =  2.54  cm. 
%  Indicates  run-time-constrained  result. 
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Fig.  8a  shows  how  the  velocity  of  the  disks  varies  with  distance.  This  plot  shows 
that  the  need  for  relativistic  equations  of  motion  is  quite  real.  Fig.  8b  is  a  blow-up  of 
Fig.  8a  which  shows  that  the  velocity  reaches  2.98x1 010  cm/s  by  the  time  a  disk  gets 
only  5  %  of  the  way  across  the  gap. 

Fig.  9a  shows  how  the  density,  similar  to  the  velocity,  is  nearly  constant  across 
the  gap  (at  equilibrium).  This  is  to  be  expected  at  equilibrium,  as  all  of  the  disks  have 
the  same  charge  density  and  most  are  moving  at  essentially  the  same  speed.  Fig.  9b 
is  a  blow-up  of  Fig.  9a,  showing  that  the  value  of  the  density  at  5  %  of  the  way  across 
the  gap  is  within  approximately  0.5  %  of  the  value  at  the  anode.  The  density  is 
greater  near  the  cathode  because  the  slower  moving  electrons  there  are  still 
"bunched",  whereas  those  farther  away  are  fully  separated  by  the  acceleration 
process. 

If  the  density  is  essentially  constant,  the  electric  field  should  appear  as  a  sloping 
straight  line  centered  about  E^p.  Fig.  10a  shows  the  magnitude  of  the  electric  field 
vs.  distance.  As  predicted,  the  plot  is  straight  and  centered  about  E^.  In  order  to 
see  the  effect  of  the  higher  density  near  the  cathode,  the  plot  was  blown-up  into 
Fig.  10b.   Indeed,  a  small  deviation  can  be  seen  near  the  cathode,  as  expected. 

Running  the  program  over  and  over  again  with  the  same  parameters  just  to  get 
plots  of  different  variables  can  be  time  consuming.  So  three  of  the  plots  were 
combined  into  one  multi-plot.  The  multi-plot  program  was  then  modified  to  make  it  into 
a  "cartoon,"  that  is,  the  plots  would  update  on  the  screen  with  each  chosen  (by  the 
programmer)  time  step.  Fig.  11a  shows  a  multi-plot  cartoon  freeze-frame  at  a  time  just 
before  the  current  density  reached  the  absolute  minimum.  As  expected,  the  density 
plot  is  high  near  the  anode  and  low  near  the  cathode.   The  electric  field  plot  shows 
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that  the  E-field  at  the  half-way  point  across  the  gap  is  less  than  E^,  also  indicating 
that  there  truly  is  more  charge  in  front  of  the  mid-point  than  behind  it. 

In  Fig.  11b  the  time  is  frozen  at  a  point  where  the  current  density  has  just 
reached  a  local  maximum.  The  density  graph  shows  high  density  at  either  end  of  the 
diode  (newly  created  at  the  cathode  and  what's  left  of  the  charge  at  the  anode  from 
Fig.  11a).  The  electric  field  curve  "kinks"  in  tandem  with  the  new  charge  in  the  density 
plot  below  it. 

Fig.  11c  shows  the  frame  in  which  the  current  density  approaches  a  local 
minimum.  Note  how  the  density  is  more  spread  out  than  in  Fig.  11a,  and  how  there 
is  less  curvature  to  the  E-field  curve.  By  Fig.  11d,  the  next  local  maximum  of  the 
current  density  has  been  reached,  and  by  now  the  density  is  nearly  equal  everywhere. 
The  electric  field  plot  is  near  linearity. 

In  Fig.  11e  the  multi-plot  is  run  to  the  steady-state.  The  density  is  constant  at  all 
points  except  near  the  cathode;  the  electric  field  curve  is  straight  and  centered  about 
Eapp,  as  expected.  Based  on  proper  behavior  of  the  above  plots,  it  can  be  concluded 
that  the  simulation  program  is  producing  plausible  data. 
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Figure  8a.  Space  charge  velocity  vs.  distance  ratio  at  equilibrium. 
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Figure  8b.    Space  charge  velocity  vs.  distance  ratio  at  equilibrium;  first  five  percent  of 
gap. 
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Figure  9a.  Space  charge  density  vs.  distance  ratio  at  equilibrium. 
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Figure  9b.    Space  charge  density  vs.  distance  ratio  at  equilibrium;  first  five  percent  of 
gap. 
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Figure  10a.  Electric  field  vs.  distance  ratio  at  equilibrium. 
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Figure  10b.  Electric  field  vs.  distance  ratio  at  equilibrium;  first  five  percent  of  gap. 
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V.   CONCLUSIONS 

As  expected,  the  final  equilibrium  value  of  the  current  density  for  a  pure  field 
emitter  was  found  to  be  less  than  the  Child-Langmuir  space  charge  limited  current 
density  for  the  non-run-time-constrained  simulations.  A  comparison  of  the  values  is 
presented  in  Table  II.  The  run-time-constrained  results  (denoted  by  a  $)  cannot  be 
used  to  either  verify  or  contradict  Eq.  (4),  but  based  on  the  arguments  of  section  IID.  it 
is  likely  that  verification  would  be  achieved  if  longer  run-times  were  feasible. 

Based  on  the  simulation  results  shown  in  Table  I,  it  can  be  stated  that,  except  in 
the  cases  denoted  by  a  J,  that  JF  <  JE.  Thus,  for  a  1  in  gap,  a  purely  field-emitted 
current  density  is  not  sufficient  to  cause  the  explosion  mechanism  of  a  whisker  on  the 
cathode  to  trigger  in  less  than  10  ns.  The  above  noted  exceptions  may  or  may  not 
conform  to  the  other  data;  computer  run  time  constraints  did  not  allow  for  a  definitive 
determination. 

Again,  it  must  be  pointed  out  that  this  model  of  only  pure  field  emission  is  quite 
elementary  in  that  it  does  not  include  the  effects  of  thermionic  emission.  A  whisker 
certainly  heats  up  due  to  the  emission  process  and  some  thermionic  emission  must  be 
present.  In  fact,  based  on  the  results  of  those  simulations  which  ran  until  equilibrium 
was  reached  (mEapp  <  3.94x1 07  V/cm),  it  would  seem  that  thermionic  emission  must 
be  the  dominant  emission  source,  at  least  for  those  parameters,  if  the  explosive  model 
is  truly  accurate.  Whether  enough  thermionic  emission  occurs  to  validate  the 
explosive  model  is  debatable,  especially  in  cases  such  as 
m  =  1,  Egpp  <  3.94x1 06  V/cm,  where  the  field  emission  contributes  almost  nothing  to 
the  total  emission  current  density  and  the  applied  electric  field  does  not  provide  much 
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relative  enhancement  (compared  to  higher  values  of  Eapp)  at  the  metal-vacuum 
interface.  A  study  done  by  Hallal  [Ref.  3]  at  such  parameters  concluded  that  the 
explosive  model  could  not  account  for  the  breakdown  that  he  observed,  and  that  a 
Schwirzke-type  model  seemed  more  promising.  It  can  be  concluded,  then,  that  the 
explosive  model  only  has  a  chance  of  validity  in  the  case  of  mE^  £  3.94x1 07  V/cm, 
with  some  small  dependence  on  work  function,  for  a  diode  with  a  1  in  gap. 
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APPENDIX  A  :  DIMENSIONLESS  VARIABLES 

Although  modern  computers  can  easily  handle  both  very  large  and  very  small 
numbers,  it  is  still  considered  good  programming  practice  to  work  with  small  ( =  1 ) 
numbers  whenever  possible.  To  facilitate  this,  variables  (such  as  time,  distance, 
velocity,  etc.)  are  transformed  into  dimensionless  quantities.  The  transformation  is 
accomplished  by  dividing  each  occurrence  of  a  particular  variable  by  some  constant 
value  of  the  variable  relevant  to  the  problem.  For  example,  if  one  were  to  examine  the 
motion  of  billiard  balls  on  a  table  top,  one  could  divide  the  masses  in  the  equations  of 
motion  by  the  mass  of  a  single  billiard  ball  mb,  thus  reducing  the  values  of  the  masses 
to  1  during  the  computations.  After  all  of  the  calculations  are  finished,  one  would  then 
multiply  the  dimensionless  values  by  the  appropriate  constants  to  get  real  values. 

Since  some  variables  are  composites  of  others  (e.g.,  velocity  is  composed  of 
distance  and  time),  once  a  certain  number  of  variables  has  been  transformed, 
transformation  of  the  remaining  ones  is  then  governed  by  the  constants  already 
chosen  and  the  physical  equations  of  the  problem  at  hand.  For  the  billiard  ball 
example,  one  might  choose  time,  distance,  and  mass  for  transformation  first. 
Momentum,  force,  and  energy,  since  they  are  composed  of  the  three  already 
transformed  variables,  would  eventually  be  transformed  using  the  three  constants  for 
time,  distance,  and  mass,  along  with  any  physical  constants  which  may  appear  in  the 
governing  equations.  Often  the  "base"  (non-composite)  variables  are  chosen  for 
transformation  first.  In  this  thesis,  such  base  variables  would  be  time,  mass,  charge, 
and  distance.  Indeed,  the  first  three  have  been  chosen  here.  But  rather  than 
distance,  current  density  has  been  selected  instead.  Current  density  was  chosen 
since  the  graph  of  current  density  vs.  time  was  so  central  in  the  thesis. 
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The  basic  transformations  are,  then: 


time  :      t  x  —  ->  x   ; 

mass  M  x  — —  ->  M  ; 

Me 

charge  :      q  x  —  ->  q   ; 


current 
density 


e 

J_ 
Jo 


All  other  variables  can  be  transformed  using  the  four  constants  above: 

2£q  Me 

distance  :       x  x  — ->  x   ; 

T3eJ0 

2  Cq  Me 
velocity  :       v  x  — - -»  v  ; 


acceleration  a  x  — : — —  -»  a  ; 


2£oMe 

T3eJ0 

2eoMe 

T2eJ0 

2eoMe 

TeJ0 

2eo 

TJf 


electric 

field  :       E  x 

Where 

T  =  time  step  (s)   ; 

Me  =  electron  mass  =  9.1 1x1 0"31  kg   ; 

e  ■  I  electron  charge  I  =  1.60x10-19  C  ; 

J0  s  current  density  at  cathode  at  t  =  0  ( A/cm* )  ;  and 

Eq  a  permitivity  of  free  space  =  8.85x10-12  C2/N-m2 


.2 
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The  extra  factor  of  2eo  comes  from  the  governing  physical  equations.  This  can 
be  seen  by  looking  at  the  derivation  of  the  transformation  of  the  dimensionless  electric 
field,  E.   Eq.  (13)  can  be  rewritten  as 

Jn  T  Jk     At  Jn  T 

E  =  T^T  =  T^JkM-  (23> 


Rearranging  gives 

2eo 


E  =  ZJkAi.  (24) 


J0T 
As  the  right  side  is  dimensionless,  so  must  be  the  left  side,  and  one  gets 

E  =  -r-TE  =  £3kAT.  (25) 

Jo  i 

If  At  =  1 ,  as  it  does  in  the  program, 

E  =  £Jk  •  (26) 

The  power  of  Eq.  (26)  lies  in  its  simplicity;  dimensionless  electric  fields  and  current 
densities  become  "equivalent."  Only  addition  is  required  for  the  calculation  of  the 
dimensionless  electric  field;  no  program  steps  are  wasted  multiplying  the 
dimensionless  current  densities  by  constants.  As  equations  of  the  form  of  Eq.  (26) 
make  up  a  large  part  of  the  computations,  the  reduction  in  the  number  of  calculations 
can  be  significant. 
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Since  the  factor  of  2eq  appears  in  the  expression  for  the  dimension  less  electric 
field,  it  must  also  appear  in  the  remaining  equations  for  acceleration,  velocity,  and 
distance.   From  Eq.  (22)  one  can  get  for  the  acceleration  a  =  dv/dt 


a  = 


j_^  E(2ep)    Me 


Y3   e      J0T 


M 


e  JnT 


2eoMe 


(27) 


Rearranging  yields 


2eoMe  gE 

eJ0T        "     y*fo 

Again,  as  the  right  side  is  dimensionless,  so  must  be  the  left  side,  and  one  gets 

a  =   lfo_Mea  =   _CL^ 
e  J0  T  y3  M 


(28) 


(29) 


Interestingly  enough,  since  q  =  M  =  1  in  the  program,  acceleration,  almost  like  current 
density,  is  nearly  "equivalent"  to  the  electric  field  in  dimensionless  form. 

Velocity  and  distance  are  related  to  acceleration  by  factors  of  (1/T)  and  (1/T)2, 
respectively: 


2  £q  Me  2  en  Me 


eJnT2 


eJnT3 


(30) 
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APPENDIX  B  :  SIMULATION  PROGRAM 


/*  Variable  Parameters  */ 


#  i  nc I ude  <s td  i  o . h> 
#include  <cmath.h> 
#include   <greeks.c> 

ne   RUN  100000 

ne   dt  1 

ne   T  1.0E-17 

ne   GAP  2.54 

ne  m  100 

ne  wkfun  4.0 

ne   Eapp  3.94E6 

ne   yscale  1000 

ne  MAX  1.1 

n  e  mu  I  t  1 

ne   hunch  3.0 

ne   I  imi  t  800 


#def 
#def 
#def 
#def 
#def 
#def 
#def 
#def 
#def 
#def 
#def 
#def 


ar  ray  s  i  ze  * / 

time  step  (d imens i on  I  ess)  /* 

t  ime  step  (sec)  * / 

gap  size  (cm)  */ 

enhancement  factor  */ 

wo  rk  f unct  i  on  (eV)  * / 

appl  ied  E-fie  I  d  (V/cm)  */ 

for  p I o 1 1  i  ng  * / 

plot  all  pts.  w/y-val.>MAX  * 

plot  limit /mu It  pts .  *  / 

ar  t  i  fie  i  al    charge   es  t  .    *  / 

#   of    run   steps    */ 
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/*  Declarations  and  Assignment  of  Initial  Values  */ 


mai  n( ) 

{ 

double   E[RUN+1 ] ,J[RUN+1 ] ,x[RUN+1 ] ,v[RUN+1 ] ,a[RUN+1 ] ,F, 
Jmu  I t  ,gapd  i  s t ,dx , xmi  n , xmax ,ymi  n ,  ymax , I abl  xmi  n , 
lab  I xmax , I abl ymi  n , I abl ymax , ploty [RUN+1 ] ,Ffact , 
J  fact  .Jexpfact  ,c2  ,  a  I  pha  ,  garrma  i  nv  ,  spchrg  ; 

int   c,f,i,j,k,l,p,q,t,o Ides t , n ,pmu I t , i  xmi  n , i  xmax , 
i  ymi  n  ,  i  ymax  ,  pi  otx  [RUN  +1  ]  .mark  ,  i  y  ,  flag  ,  i  I  y  ; 


pmu  It  =  mu  I  t  ; 

Jexpfact  =  ( -6.42E7)*pow(wkf un, 1 .5) ; 

F  =  m*Eapp ; 

J[0]  =  ((1 .54E-6)*F*F/wkfun)*exp(Jexpfact/F) ; 

Ffact  =  m*T*J[0]/(2.0*(8.85E-14)) ; 
Jfact  =  (1  .54E-6)/(wkfurTJ[0])  ; 

E[0]  =  Eapp*2.0*(8.85E-14)/(T*J[0]) ; 

alpha  =  (1 .0E-6)*2.0*(8.85E-12)*(9.11E-31 )/ 
(T*T*T*(1 .6E-19)*J[0]) ; 

gapdist  =  alpha*GAP; 

x[0]  =  0.0; 
v[0]  =  0.0; 

c2  =  alpha*T(3.0E10)*alpha*T*  (3.0E10)  ; 

o Ides t  =  1  ; 

t  -  1  ; 

mark  =  1 ; 

c  =  1  ; 

flag  =  0; 

f  =  0; 

spchrg  =  hunch; 

ymi  n  =  -0 .001 ; 
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/*  Determine  Electric  Fields,  Accelerations,  Velocities,  and 
Positions  for  Existing  Disks;  Eliminate  Disks  at  Anode  */ 


wh  i  I  e( t  <=  limit) 

{ 

for(i    =   oldest;    i    <=    (mark-1);    i    =    i    +1 ) 

{ 


i  f  ( i     !=   oldest) 

E[i]    =   E[i-1]    -    J[i-1]    -    J[i]; 


garrmainv   =   sqr  t  ( 1  .  0- ( v  [  i  ]  *  v  [  i  ] /c2) )  ; 


a[ 
v[ 
x[ 


]    =   E  [  i  ]  *  g  anrma  i  n  v  *  g  arrma  i  n  v  *  g  arrma  i  n  v  ; 
]    =   v[i]    +   a[ i ] *dt ; 
]    =   x[i]    +   v[i]*dt; 


if(x[i]    >=   gapd  i  st ) 

{ 

o I des  t    =   o I des t    +    1 ; 

for(j  =  oldest;  j  <=  (mark-1);  j 
E[j]  -  E[J]  +  J[oldest  -  1  ]  ; 

} 


=  j+1) 
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/*  Determine  Electric  Fields,  Velocity,  Position,  and 
Current  Densities  for  Newly  Created  Disks  at  Cathode; 
Decrease  Artificial  Space  Charge  */ 


i  f  (flag==0) 

E[mark]  =  E[0]  -  spchrg; 


e  I  se 


E[mark]  =  E[0] ; 


x [mark]  =  x [0] ; 
v [mark]  =  v [0] ; 


for(k  =  oldest;  k  <=  (mark-1);  k=  k+1 ) 
E[mark]  =  E[mark]  -  J[k] ; 


F  =  Ffact*E[mark] ; 

J[mark]  =  J f act *F*F*exp( Jexpf ac t /F) ; 

spch  rg  =  spch  rg  -  J  [  1  ]  ; 

i  f  (spchrg<0  .0   &&   flag   ==   0) 

{ 

flag   =    1  ; 
f    =    t; 

} 
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/*  Determine  Points  for  Plotting;  Calculate  Electric 
Field  Changes  Due  to  Newly  Created  Disk  at  Cathode; 
Increment  Time  Step  */ 


i  f ( t==pmul t ) 

{ 

pi otx [c]  =  t  ; 

plotyjc]  =  I og1 0 ( J [mark] )  ; 

c  =  c+1  ; 

pmu  I t  =  pmu  It  +  mu  I t  ; 

} 
else  i f (J[mark]>MAX) 

{ 

pi otx [c]  =  t  ; 

plotyjc]  =  I og1 0 ( J [mark] )  ; 

c  =  c  +  1  ; 

) 


i  f ( I og1 0( J [mar k] )<ymi  n ) 

ymi  n  =  I og1 0( J [mark] ) ; 


f o r ( I  =  oldest;  I  <=  (mark-1);  I  =  1+1) 
E[ I ]  =  E[ I ]  +  J [mark] ; 


t  =  t  +  dt  ; 

ma  r  k  =  ma  r  k  +  d  t  ; 
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/*  "Recycle"  Arrays  if  Full;  End  Wh i  I e  Loop  */ 

if(mark>RUN  &&  RUN !  =  l  imi  t ) 

{ 

for  (q=o  Iciest  ;q<=RUN;q=q+1 ) 

{ 

E[q-oldest+1 ]  =  E[q] ; 
J [q-o ldest  +  1 ]  =  J [q]  ; 

a[q-o ldest+1 ]  =  a[q] ; 
v [q-o ldest  +  1 ]  =  v[q]  ; 
x [q-o ldest  +  1 ]  =  x [q]  ; 

} 

mark  =  RUN  -  oldest  +2; 
o Ides  t  =  1  ; 

} 
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/ *  Setup  Graph  i  cs  * / 


ymax  =  0.0; 

i  xmi  n  =  0 ; 

i  xmax  =  limit; 

iymin  =  (int)(ymin*yscale); 
i ymax  =  ( i n t ) (ymax*yscal e) ; 

I abl xmi  n  =  0.0; 

I abl xmax  =  T* I imi  t ; 

I abl ymi  n  =  ymi  n ; 
I abl ymax  =  0.0; 

************************************************* 

ini  t()  ; 

color_scale("cyanblue") ; 
grey_scale("greyscale1" );  wi  ndow0( ) ; 
bgcol(7);  erase();  color(O); 


move 
move 
s  i  ze 
move 

move 
move 
move 
move 

move 
move 
move 
move 

move 
move 


245,30) ;pr  int f 

75,384) ;pr  intf 

2,2); 

435  ,700) ; p r  int 

85,96) ;pr  i  nt  f ( 
85,660) ;pr  int  f 
195  ,682) ; p r  int 
760  ,682) ; p r  int 

810,150) ;pr  int 
810,250) ;greek 
810,350) ;greek 
810,450) ;pr  int 


810  ,550) ;pr  int 
810,630) ;pr  int 


move(350,760) ;pr  int 
print 
sup(  ' 


"Log  Current  Density  Ratio  vs.  Time"); 
"Log  J/J")  ;size(1 ,1 )  ;sub( '0' ) ; 

("Time  (sec) ")  ; 


%-3.2f " , I abl ymax)  ; 
"%-3.2f " , lab Iymin) 
("%3.2f " , lablxmin) 
("%3.2e" , I abl xmax) 


("m  =  %d"  ,m)  ; 

BETA) ;pr intf ("=%-2.1e" ,T) ; 

Phi  )  ;pr  int  f ("  =  %-2  . 1 f " ,wkf un) ; 

("E") ;size(1 ,1 ) ;sub( "a* ) ;sub( 'p* ) ; 


sub( ' p" ) ;size(2,2) ;pr  int f ("=%-2. 1e" ,Eapp) ; 


("f  =  %d",f)  ; 
("old=%d" .oldest) ; 

("J")  ;size(1 ,1) ;sub('0* ) ;size(2,2) ; 
("  =  %4.2e  Amp-cm" ,J[0]) ;size(1 ,1 ) ; 
,);sup(,2,);size(2,2); 


wi ndow(224, 96 ,800,672) ;bg col (3) ;e rase O ; 

sea  I e ( i  xmi  n ,  i  ymax , i  xmax , i  ymi  n ) ; 

rect ( i  xmi  n , i  ymax , i  xmax , iymin) ;color(4) ; 
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/*  Plot  Graph  of  Log  Current  Density  vs  Time;  End  Program  */ 


catvec( ixmin, ( int) ( I og 1 0 ( J [ 1 ] ) *ysca  le) )  ; 

for (p=1 ;p<c;p=p+1 ) 

pr int f ("%d,%d,0, plot x[p] , ( int) (pi oty [p] 'yscal e) ) ; 

P  r  i  n  t  f  ( "  ; " )  ; 


windowO( );color(0);size(1,1) ;mode( "PO" ) ; 

iy  =  672-(int)((672-96)*((yscale*ploty[p-1])-iymin)/ 
( i  ymax - i  ymi  n) ) ; 

move(805  ,  i  y )  ; 

p  r  i  n  t  f ( "%4  .  3e " , p I o  t  y [ p - 1 ] ) ; 

ily  =  672-( int)((672-96)*((yscale*ploty[1])-iymin)/ 
( i  ymax - i  ymi  n) ) ; 

move( 1 60  ,  i  I y )  ; 

pr intf ("%4.3e" ,ploty[1 ] ) ; 

} 
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